#include <stdio.h>
#include <iostream>
#include <fstream>
#include "matrix2.h"


//------------------------------------------------//
//------------------CONSTRUCTORS------------------//
//------------------------------------------------//

//----------------------------------------------------------------//
//Constructor
//@param vals : values of the scoring values used
//----------------------------------------------------------------//
Matrix2::Matrix2(float * vals)
{
	match_score=vals[0];
	mismatch_score=vals[1];
	gSeq_score=vals[2];
	gMono_score=vals[3];
	bestA_score=-1;
	score=NULL;
	dir=NULL;
	match=NULL;
	align_start_x=align_end_x=align_start_y=align_end_y=-1;
	align_len=0;

}




//------------------------------------------------//
//-------------------DESTRUCTOR-------------------//
//------------------------------------------------//
Matrix2::~Matrix2()
{
	int i;
	if(score!=NULL)
	{
		for(i=0; i<size_x; i++)
		{
			delete score[i];
			delete dir[i];
			delete match[i];
		}
		delete [] score;
		delete [] dir;
		delete [] match;
	}
}




//------------------------------------------------//
//------------------FUNCTIONALITY-----------------//
//------------------------------------------------//

//----------------------------------------------------------------//
//Create and populate matrices and call traceback()
//@param monoSeq : the monomer sequence
//@param monoSeqLen : length of the monomer sequence
//@param subSeq : subsequence from the search sequence
//@param subSeqLen : subsequence from the search sequence length
//@param o : start location of  subsequence within search sequence
//@return bool : true
//----------------------------------------------------------------//
bool Matrix2::make_matrices(char *&monoSeq, int monoSeqLen, char *&subSeq, int subSeqLen, int o)
{
	int i,j;

	size_x=monoSeqLen;						//monomer always x
	size_y=subSeqLen;
	offset=o;

	//this will store the best alignment for the matrix
	bestA_score=0;

	//set up arrays
	dir=new char * [size_x];
	score=new float * [size_x];
	match=new bool * [size_x];
	for(i=0; i<size_x; i++)
	{
		dir[i]=new char [size_y];
		score[i]=new float [size_y];
		match[i]=new bool [size_y];
	}
	
	//fill matrices
	for(i=0; i<size_x; i++)
	{
		for(j=0; j<size_y; j++)
		{
			if(monoSeq[i]==subSeq[j])
			{
				match[i][j]=true;
				score[i][j]=calc_score(i,j,dir[i][j]);
			}
			else
			{
				match[i][j]=false;
				score[i][j]=calc_score(i,j,dir[i][j]);
			}
		}
	}
	traceback();
	return true;
}

//----------------------------------------------------------------//
//Determines best alignment from matrix values and stores the
//start and end positions of alignment found in class vars
//@return bool : true
//----------------------------------------------------------------//
bool Matrix2::traceback()
{
	//traceback
	int i,j,s;
	float max_val=0;
	int number_of_start_points;
	bool loop;

	max_val=get_max_value();
	number_of_start_points=get_num_starting_pos(max_val);


	//int align_start_x, align_end_x, align_start_y, align_end_y;
	//these are temporary values for the traceback... 
	//in the even that there is more than 1 starting point
	//temporary values are needed.
	int t_align_end_x=0;
	int t_align_end_y=0;
	int t_align_start_x=0;
	int t_align_start_y=0;

	float align_score;
	int t_align_len;

	for(s=0; s<number_of_start_points; s++)
	{
		//determine index of max_value & set the first node of the path
		for(i=0; i<size_x; i++)
		{
			for(j=0; j<size_y; j++)
			{
				if(score[i][j]==max_val)
				{
					t_align_end_x=i;
					t_align_end_y=j;
					score[i][j]=0;			//makes sure I don't pick this value as a max again
					i=size_x;				//gets me out of the loop to find the start point
					j=size_y;
				}
			}
		}

		//wander back from this point... accumulate score as you go
		t_align_len=0;
		align_score=0;
		loop=true;
		i=t_align_end_x;
		j=t_align_end_y;
		while(loop)
		{
			if(match[i][j])
				align_score++;
			switch(dir[i][j]){
			case 'U': i--; break;
			case 'L': j--; break;
			case 'S': i--; break;
			default: i--; j--; break;		//D, Q, R, or T
			}
			t_align_len++;

			if( ((i==0) && (j==0)) || (i==0) )
				loop=false;
		}
		
		align_len=t_align_len+1;
		align_score=(align_score/(float)t_align_len)*100.0;
		
		if(align_score>bestA_score)
		{
			bestA_score=align_score;
			align_start_x=t_align_start_x+1;
			align_end_x=t_align_end_x+1;
			align_start_y=t_align_start_y+1;
			align_end_y=t_align_end_y+1;
		}
	}
	return true;
}

//----------------------------------------------------------------//
//Gets the maximum score in the score matrix
//@return float : max_val
//----------------------------------------------------------------//
float Matrix2::get_max_value()
{
	int i,j;
	float max_val=0;
	for(i=0; i<size_x; i++)
	{
		for(j=0; j<size_y; j++)
		{
			if(score[i][j]>max_val)
			{
				max_val=score[i][j];
			}
		}
	}
	return max_val;
}

//----------------------------------------------------------------//
//Determines the number of starting positions with the best
//alignments according to the score matrix
//@param max : maximum score in the score matrix
//@return int : num_start
//----------------------------------------------------------------//
int Matrix2::get_num_starting_pos(float max)
{
	int i,j;
	int num_start=0;
	for(i=0; i<size_x; i++)
	{
		for(j=0; j<size_y; j++)
		{
			if(score[i][j]==max)
			{
				num_start++;
			}
		}
	}
	return num_start;
}

//----------------------------------------------------------------//
//Determines the individual cell scores and direction for
//populating the score and direction matrix
//@param x : cell position row
//@param y : cell position column
//@param d : direction - pass by reference
//@return float : the score
//----------------------------------------------------------------//
float Matrix2::calc_score(int x, int y, char & d)
{
	float up,left,diag;

	//diagonal
	if(match[x][y])
	{
		if( (x>0) && (y>0) )
			diag=score[x-1][y-1]+match_score;
		else
			diag=1;							//when initializing to zero
	}
	else
	{
		if( (x>0) && (y>0) )
			diag=score[x-1][y-1]+mismatch_score;
		else
			diag=0;							//when initializing to zero
	}

	//up   --- gap in sequence
	if(x!=0)
		up=score[x-1][y]+gSeq_score;
	else
		up=0;

	//left  --- gap in monomer
	if(y!=0)
		left=score[x][y-1]+gMono_score; //(score[x][y-1]/(score[x][y-1]+1));
	else
		left=0;

	if(x==0)	//must be left
	{
		d='L';
		return left;
	}
	
	if(y==0)	//must be up
	{
		d='U';
		return up;
	}

	//determine maximum
	if( (up==0) && (left==0) && (diag==0) )
	{
		d='D';
		return diag;
	}
	if( (diag>up) && (diag>left) )
	{
		d='D';
		return diag;
	}
	if( (up>diag) && (up>left) )
	{
		d='U';
		return up;
	}
	if( (left>diag) && (left>up) )
	{
		d='L';
		return left;
	}
	if( (diag==up) && (diag>left) )
	{
		d='Q';				//diag or up
		return diag;
	}
	if( (diag==left) && (diag>up) )
	{
		d='R';
		return diag;		//diag or left
	}
	if( (left==up) && (left>diag) )
	{
		d='S';
		return left;		//up or left
	}
	
	d='T';					//default = diag or left or up
	return diag;
}




//------------------------------------------------//
//-------------------ACCESSORS--------------------//
//------------------------------------------------//

//----------------------------------------------------------------//
//Retrieve alignment from vals stored in this
//@param a : align instance to create
//@return bool : true if matrices have been populated
//----------------------------------------------------------------//
bool Matrix2::get_alignment(Align *& a)
{
	if(score==NULL)
		return false;
	if(a!=NULL)
		delete a;
	a=new Align(align_start_x,align_end_x,align_start_y+offset,align_end_y+offset,align_len,bestA_score);
	return true;
}


float Matrix2::get_best_align_score() {return bestA_score;}


